#include <stdio.h>
#include <stdlib.h>

typedef struct
{
    int J_ORDER;            // Ordem da matriz
    int J_ROW_TEST;         // Linha para ser avaliada
    double J_ERROR;         // Erro permitido
    int J_ITE_MAX;          // Número máximo de iterações
    double **MA;            // Matriz A do sistema         
    double *MB;             // Matriz B do sistema

    double *MA_J_ROW_TEST;  // Armazena a linha a ser avaliada
    double MB_J_ROW_TEST;   // Armazena o resultado da linha a ser avaliada
}Jacobi;


/*
*   Aloca as matrizes da struct 'Jacobi' e inicializa seus valores pela leitura do arquivo.
*/
void inicia_Jacobi(FILE *pFile, Jacobi *J)
{
    int i, j;

    fscanf (pFile, "%d", &(J->J_ORDER));
    fscanf (pFile, "%d", &(J->J_ROW_TEST));
    fscanf (pFile, "%lf", &(J->J_ERROR));
    fscanf (pFile, "%d", &(J->J_ITE_MAX));

    /* Aloca a matriz MA */
    J->MA = malloc(J->J_ORDER * sizeof(double));

    for(i=0;i<J->J_ORDER;i++)
        J->MA[i] = (double *) malloc(J->J_ORDER * sizeof(double));

    /* Aloca o vetor MB */
    J->MB = (double *) malloc(J->J_ORDER * sizeof(double));

    /* Leitura da matriz MA */
    for(i=0;i<J->J_ORDER;i++)
        for(j=0;j<J->J_ORDER;j++)
            fscanf (pFile, "%lf", &(J->MA[i][j]));

    /* Leitura do vetor MB */
    for(i=0;i<J->J_ORDER;i++)
        fscanf (pFile, "%lf", &(J->MB[i]));

    J->MA_J_ROW_TEST = (double *) malloc(J->J_ORDER * sizeof(double));
}

/*
*   Prepara as matrizes MA e MB para a aplicação do método de Jacobi-Richardson.
*
*   A linha de teste(J_ROW_TEST) é armazenada na struct para o teste final de checagem.
*
*   Cada linha das matrizes, então, é dividida pelos respectivos valores da diagonal da mesma linha.  
*
*   No vetor 'x' é armazenado os valores iniciais do método Jacobi-Richardson. 
*/
void prepara_Jacobi(Jacobi *J, double *x)
{
    int i, j;
    double div;

    for(i=0;i<J->J_ORDER;i++)
    {
        /* Armazena a linha caso seja a linha do teste final (J_ROW_TEST) */
        if(i == J->J_ROW_TEST)
        {
            for(j=0;j<J->J_ORDER;j++)
            {
                J->MA_J_ROW_TEST[j] = J->MA[i][j];
            }

            J->MB_J_ROW_TEST = J->MB[i];
        }

        /* Armazena a diagonal da linha */
        div = J->MA[i][i];

        /* Divisão da linha pela diagonal */
        for(j=0;j<J->J_ORDER;j++)
        {
            J->MA[i][j] /= div;
        }
        J->MB[i] /= div;

        /* Zera a diagonal da linha pois é ignorada no cálculo */
        J->MA[i][i] = 0;

        /* Armazena no vetor 'x' os valores iniciais do método Jacobi-Richardson */
        x[i] = J->MB[i];
    }
}

/*
*   @Método Jacobi-Richardson
*
*   'x' é uma matriz de 2 linhas e J_ORDER colunas. 
*
*   A partir dos valores anteriores armazenados em x[id][...], realizamos uma nova iteração
*   e armazenamos os novos valores novamente em x[1][...] ou x[0][...], dependendo do valor de 'id', que sempre vale 0 ou 1.
*
*   O erro a partir dos valores de 'x' é então retornado.  
*/
double jacobiRichardson(Jacobi J, double **x, int *id)
{
    int i, j;
    double aux;
    double erro1 = 0;   // Numerador no cálculo do erro final
    double erro2 = 0;   // Denominador no cálculo do erro final

    /* Os valores novos devem ser armazenados x[1][J_ORDER] */
    if(*id == 0)
    {
        for(i=0;i<J.J_ORDER;i++)
        {
            x[1][i] = J.MB[i];

            /* Realiza as devidas multiplicações e subtrações da linha */
            for(j=0;j<J.J_ORDER;j++)
            {
                x[1][i] -= x[0][j] * J.MA[i][j];
            }

            /* Armazena o erro */
            aux = x[1][i] - x[0][i];

            if(aux < 0) aux *= -1;

            if(aux > erro1 ) erro1 = aux;

            aux = x[1][i];

            if(aux < 0) aux *= -1;

            if(aux > erro2 ) erro2 = aux;
        }

        *id = 1;
    }
    /* Os valores novos devem ser armazenados x[0][J_ORDER] */
    else
    {
        for(i=0;i<J.J_ORDER;i++)
        {
            x[0][i] = J.MB[i];

            /* Realiza as devidas multiplicações e subtrações da linha */
            for(j=0;j<J.J_ORDER;j++)
            {
                x[0][i] -= x[1][j] * J.MA[i][j];
            }

            /* Armazena o erro */
            aux = x[0][i] - x[1][i];

            if(aux < 0) aux *= -1;

            if(aux > erro1 ) erro1 = aux;

            aux = x[0][i];

            if(aux < 0) aux *= -1;

            if(aux > erro2 ) erro2 = aux;
        }

        *id = 0;
    }

    /* Retorna o erro */
    return erro1/erro2;
}

/*
*   Imprime a saída pelo padrão do projeto.
*/
void saida(Jacobi J, double **x, int id, int iter)
{
    int i, j;
    double aux=0;

    for(i=0;i<J.J_ORDER;i++)
    {
        aux += J.MA_J_ROW_TEST[i] * x[id][i];
    }

    printf("Iterations = %d\n", iter);
    printf("RowTest: %d => [%lf] =? %lf", J.J_ROW_TEST, aux, J.MB_J_ROW_TEST);
}

int main()
{
    int i;
    int iter = 0;   // Número de iterações
    int id = 0;     // Índice para controlar os resultados anteriores e novos da matriz 'x' (Será sempre 0 ou 1) ´
    double **x;     // 'x' irá guardar os novos valores e os anteriores de cada iteração do método numérico 
    double error;   // Armazena o erro do método 
    FILE *pFile;
    Jacobi J;

    pFile = fopen ("matriz500.txt","r");

    /* Inicializa a struct */
    inicia_Jacobi(pFile, &J);

    /* Aloca x onde será armazenado os resultados das iterações */
    x = malloc(2 * sizeof(double));

    for(i=0;i<2;i++)
        x[i] = (double *) malloc(J.J_ORDER * sizeof(double));

    /* Prepara a struct */
    prepara_Jacobi(&J, x[id]);

    /* @Jacobi-Richardson */
    do
    {
        error = jacobiRichardson(J, x, &id);

        iter++;

    /* Sai do método caso o erro seja inferior a J.J_ERROR ou ultrapasse o número máximo de iterações permitidas*/
    }while(error > J.J_ERROR && iter < J.J_ITE_MAX);

    /* Imprime a saída padrão */
    saida(J, x, id, iter);

    printf("\n\nFIM");

    return 0;
}
